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ABSTRACT 


The bulge in the Earth at its equator has been shown [1] 
to lead to a clustering of natural decays biased to occur 
towards the equator and away from the orbit’s extreme 
latitudes. Such clustering must be considered when 
predicting the Expectation of Casualty (E.) during a 
natural decay because of the clustering of the human 
population in the same lower latitudes. This study 
expands upon prior work [1, 2], and formalizes the 
correction that must be made to the calculation of the 
average exposed population density as a result of this 
effect. Although a generic equation can be derived from 
this work to approximate the effects of gravitational and 
atmospheric perturbations on a final decay, such an 
equation averages certain important subtleties in 
achieving a best fit over all conditions. The authors 
recommend that direct simulation be used to calculate the 
true E, for any specific entry as a more accurate method. 
A generic equation is provided, represented as a function 
of ballistic number and inclination of the entering 
spacecraft over the credible range of ballistic numbers. 


1. BACKGROUND 


In the final stages of orbital decay, the density scale 
height shrinks to only a few kilometers. At the same 
time, the rapidly-circularizing orbit encounters two 
effects of the oblate Earth, each of which (at moderate to 
high inclinations) permute the “effective altitude” to a 
degree comparable to or exceeding such typical scale 
heights. The primary influence is the elevation of the 
atmosphere relative to Earth’s center in the equatorial 
regions, simply rising in lock step with the geoid. 
Because the Earth’s radius at the equator is 21 kilometers 
(km) higher than at its poles, the effective density (and 
thus, the local deceleration) could hypothetically increase 
by a factor of nearly 50 at the equator relative to values 
at the poles under an idealized perfect circular Kepler 
orbit. 


Simultaneously, the J2 gravitational orbit perturbation 
works in the opposite direction, raising the local radius of 
the orbit at extreme latitudes, and dropping it at the 
equator. There are additional J3 effects noted by 
Fremeaux, et al., [3] which emerge as a small correction 
factor in the present study. Fig. 1 shows the progression 
of the local altitude and radius for a decaying polar orbit, 
and the progression in the Logio of density during the 


decay. All values are plotted versus latitude, causing a 
rocking back and forth as time progresses. Clearly, the 
dominant nonlinear factor in the otherwise spiral decay is 
the physical uplift of the atmosphere into the orbit path. 


Minus Log Density 
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Figure 1. Altitude, radius, and Minus Log(Density) plotted 
vs. latitude for a decaying 50° orbit at 100kg/m2 ballistic 
number. While the J2 gravitational effect raises the local 

radius at the extreme latitudes, the elevation of the 
atmosphere above the physical bulge at the equator 
dominates the density variations encountered by the 
decaying spacecraft. Both the gravitational and atmospheric 
perturbations are evident in structure of the decay 
concentration curves. 


With such strong cyclic variation in the drag profile, it 
should be expected that the distribution of natural decays 
around the orbit path (measured in Argument of Latitude 
[ArgLat with symbol ©], which is the angle along an orbit 
in the orbit plane measured from the north-bound pass of 
the equator; generally computed as the true anomaly plus 
the argument of perigee) should instead show a statistical 
clustering that reflects this cyclic forcing function at one- 
half orbit period. Indeed, in earlier work, both simulation 
and analysis of well-documented real, natural polar orbit 
decays showed a cyclic clustering towards the equator. 
Preliminary assessments with the initial model revealed 
that the statistical biasing of natural decays to the more 
heavily populated lower latitudes could have up to a 13% 
effect on the calculated value of E.. When the authors 
presented the preliminary model they acknowledged that 
further work would be necessary to completely define the 
Statistical biases under different values of ballistic 
number (BN), inclination, and perhaps solar flux 
conditions. This paper reports the conclusion of that 
work, covering a full range of credible ballistic numbers 


and all orbit inclinations. Solar flux was shown in 
subsequent publication [2] to be inconsequential to the 
decay latitude bias problem. In the improved model the 
authors have reformulated the analysis in terms of ArgLat 
instead of latitude, which has allowed better comparison 
of the physics across cases. 


2. APPROACH 


Tens of thousands of natural decays were simulated in the 
NASA General Mission Analysis Tool (GMAT), in 
parametric studies every 5° of inclination for spacecraft 
with 50, 100, 150, and 200 kg/m? BN decaying for 
slightly more than 50 orbits each. The orbits were 
propagated to 90 km altitude. (Past research [2] has 
shown that at 90 km the residual time and distance to go 
to the surface is approximately independent of ArgLat for 
any given ballistic number. Therefore, to reduce 
computation time from weeks to days, the decays were 
truncated at 90 km). This truncation must be accounted 
for when considering actual expectation of casualty, 
because ground intercept will be several degrees 
downrange of the calculated entry point. 


GMAT is a certified NASA mission planning tool with 
an adaptive step, ninth order Runge-Kutta integrator with 
eighth order error control that can numerically integrate 
the path of a decaying object through a selectable 
atmosphere. In this study, the MSIS-E 2000 atmosphere 
was modeled at F(10.7) identical daily and 90-day 
average values of 115 janskys, and the global 
geomagnetic index Kp constant was uniform at 3. Orbits 
were propagated with 4th-order spherical harmonics and 
full atmospheric responses to day-of-year and the orbit’s 
Right Ascension of the Ascending Node (RAAN). 


The mass of the modeled spacecraft was dithered in tiny 
steps to achieve a cumulative 2% spacecraft mass 
variation in 101 uniform steps per run. The mass dither 
is considered by the authors to be the subtlest and most 
linear way to induce a smooth variation in a very 
nonlinear process. This method allows the initial 
conditions to be as nearly identical as possible, while 
allowing the forces being studied to integrate in the 
nonlinear ways under study. The 2% total mass variation 
in a decay with more than 50 remaining orbits at 
inception is enough to observe one full cycle of the 
ArgLat of decay. (One hundred and one runs are needed 
to generate 100 gaps. It is the gaps that drive the 
compression data.). 


Once a full cycle of ArgLat was observed, the relative 
compression (or “compression factor’) C(@) of the 
decays was recorded at each value of ©. C(O) is 
calculated by talking the ratio of the spacing between the 
ArgLat of any one iteration’s decay and the next 
iteration’s value, divided by the averaged step size over 
the full 360° cycle of ArgLat. (The non-integer number 
of steps is linearly interpolated between integer steps that 


bracket the 360° cycle.) A clustering of decays in one 
portion of the arc then yields a ratio less than one, and a 
rarefaction or sparseness of decays will lead to a number 
greater than one. Although this may seem 
counterintuitive, representing a variable angle per 
incremental decay scenario (vs. variable decays per unit 
angle) was the most easily immediately-normalized 
function when establishing common compression curves. 
When the function is inverted to a probability distribution 
function, each local point must be binomially expanded, 
as discussed later. 


Example: if 57.137 uniform increments of a 
satellite’s mass caused the decay’s ArgLat to vary 
by 360°, the average step size would be 
(360°/57.137) = 6.2751°. If the argument of 
decay latitude shifted only 5.72° from one 
iteration to the next at a particular region, the local 
compression rate would be (5.72°/6.2751) = 
0.9115° 
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Figure 2. An example relative compression curve for 
BN=200 decays at 50° inclination and RAAN=180° The 
curve illustrates the excellent fit of the five-term Fourier 
approximation. At the vernal equinox, this orbit has a beta 
angle of zero and passes directly under the subsolar point, 
creating a day-night asymmetry in the decay rate. The near 
mirror of this curve occurs at 90° RAAN (see Fig 6). The 
asymmetry is removed by plotting the rates at four 
orthogonal RAANS and averaging the curves. (This is best 
done after curve-fitting, because there is no way to align the 
individual statistical points to achieve the same values of 
ArgLat in all four plots.) Note in the curve how sparse the 
entry events are at the peaks relative to the valleys, in terms 
of spatial density of the points. 


The resulting plot of relative compression vs. ArgLat 
shows a consistent cyclic behavior of concentrating 
entries as the decay location approaches the equator 
(ArgLat= 0 and 180), receding rapidly after passing it, 
and becoming most rarefied approaching the extreme 
latitudes (ArgLat = 90 and 270). An example plot is 
shown in Fig. 2. A 5-term Fourier analysis of the data is 


superimposed as the solid line, showing an excellent fit, 
with generally <0.7% average error of the fit from 
recorded values. The universal model of the compression 
curve over all conditions is generated by characterizing 
the trend of the first four best-fit Fourier coefficients 
under varying ballistic number, RAAN, and inclination. 
(The fifth term represents residual noise, but no implicit 
physics emerges from its trending.) Noise comparable to 
or greater than the fifth Fourier amplitude coefficient is 
introduced in the estimation of the curve in 256 discreet 
bins of ArgLat with nearest value in the true compression 
curve. 


The entire series was repeated at four orthogonal orbit 
planes (RAAN=0, 90, 180, and 270) at the vernal equinox. 
The reason for RAAN variation was to average out the 
effects of the diurnal variation in the atmosphere, which 
presents variations in density of up to 25% of peak value, 
depending upon whether @ lies near the terminator, subsolar 
point, or solar midnight. 


In all then, (4 ballistic numbers) x (4 RAANs) 
x (101 runs/case) x (18 inclinations) led to 29,088 sample 
runs, leading to 288 curves of 100 data points each. An 
additional 30 special cases were run to explore conditions at 
the solstice, simplified gravitational models, and extreme 
ballistic numbers, for yet another 3030 runs. 


Next, all 100 compression data pairs for each set of input 
conditions were spread into the closest of 256 uniform 
ArgLat bins. A Fourier analysis was performed on this 
256-bin representation of the curve. The amplitude and 
phase of the first five oscillating terms were recorded, with 
the first term (the uniform average compression) defined as 
being one despite any curve-fitting residual. 


The five Fourier terms were then explored for empirical 
trends, such that the amplitude and phase of each term can 
be derived for any input. The empirical formulation of the 
Fourier terms then, in theory, allows us to recreate the 
compression curve for any condition. Ultimately, chaotic 
trends in the phases of asymmetric Fourier terms, and low- 
inclination orbits in general were problematic to developing 
a fully-universal empirical prediction. 


From the relative compression C(®) along such a 
reconstituted curve we can derive the relative likelihood of 
impacting at any particular latitude. For any given 
inclination (i) and ArgLat (©), we get that the latitude is 
ARCSIN(SIN(G)*SIN(®)), and the local _ statistical 
likelihood at that latitude relative to unit average will be a 
binomial expansion of 1/C. The coefficient 1/C must be 
binomially expanded to achieve the original true 
compression in relative decays per degree. This problem 
stems from the mathematical inequality that the most easily 
normalized initial compression curve is based upon 
N/Avg(X), which is not equal to N*Avg(1/X). 


3. RESULTS 


Each compression curve takes the general form of a 
sinusoidal curve in (two times the) ArgLat, (a 
2© harmonic), perturbed by 10, 30, and 40 harmonic 
terms of generally lower amplitude, each term affected 
by a unique phase offset. 
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Figure 3. A family of compression curves for inclinations 
from 5°to 90° at RAAN=0° and BN=100 kg/n?. Curve 
amplitude near 90° ArgLat increases with orbit inclination. 
The curves are more chaotic and invert the compression 
effect’s locations at the lower inclinations, although at low 
amplitude and low inclination, the net effect on latitude 
clustering is small. 
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Figure 4a. The even-numbered (symmetric) Fourier term 
amplitudes of compression curves. Note that the 20 terms 
are typically a factor of 10 larger than the 40 terms and the 
patterns generally repeat with comparatively small 
variations for ballistic number and RAAN. In each of the 
16 tooth-shaped curves, orbit inclination increases from 5 at 
the left to 90° at the right. Each of four major groupings 
shows the results for one RAAN value: 05 905 1805 and 
270° from left to right. Within each group, each of four 
ballistic numbers are evaluated: 50, 100, 150, and 
200 kgf’. 
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Figure 4b. The odd (asymmetric) Fourier term amplitudes of 
the compression curves. The groupings are as in Fig. 4a, 
and the vertical axis is to the same scale. The simple first 

harmonic wave in @ is the largest but most irregular of all 
the secondary terms, and is believed to be heavily influenced 
by the exact path through or near the subsolar point, 
potentially containing some J3 gravitational information 
too. Note the interesting and as yet unexplained gap in most 
10 amplitudes near 50~-55° The 3 terms are very 
repeatable and are likely purely atmospheric in nature (see 
Fig 5c). The 50 terms are completely negligible, 
representing noise. (The gravitational model itself is only to 
four harmonic terms.) 
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Figure 5a. The scale of atmospheric vs. gravitational effects 
is evident in the calculation of compression curves when 
Earth’s simple point-mass gravitational model (0-order, 
0-degree gravitational model) is used instead of the full 
4-order, 4-degree gravitational model. Only atmospheric 

geometry—especially the solar heating effect—is driving the 

O, 30, and 40 Fourier terms in this case, resulting in 
perturbing effects generally comparable in scale with the 
“wall of air” effect at low inclinations, dropping to about 
25% influence at the highest inclinations for the 10 term, 
and <5% influence for the higher-order terms. 
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Figure 5b. The compression curve Fourier coefficients when 
the full gravitational model is used for the same conditions 
of 100 kg/m2 ballistic number and 0° RAAN. Note the 
reduction in nearly all of the coefficients, indicating that 
gravitational forces in many cases may mitigate the 
atmospheric effects. 
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Figure 5c. The point mass gravitational model for a 
100 kg/m? BN at 90° RAAN. Note the shift in height and 
shape of the 10 (blue) term, and the increase in the 
20 (orange) term, all presumably due to atmospheric effects 
only. The path of the orbit relative to the subsolar point has 
a noticeable effect on the symmetry of the curve (Fig. 6). 
Note that the 30 (red) term is remarkably consistent across 
all the curves. This strongly implies that the J3 terms are 
atmospheric in nature, and real. 


The basic 2© curve has an uprange offset of 15°-20° 
relative to nodal crossings at 0° and 180° of ArgLat at 
higher inclinations, with interesting trends in the opposite 
direction at the lower inclinations, where all the Fourier 
terms have comparable amplitudes. The shape of the 
curve is generally more disorganized at lower inclination 
orbits, but still exhibits this basic behavior (see Fig. 3). 


The first five Fourier terms’ amplitudes are shown in 
Figs. 4a and 4b grouped into terms that are respectively 
even and odd multiples of ©, representing the effects that 
are symmetric north-to-south and those that are 
asymmetric. 
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Figure 5d. The same conditions as in Fig. 5c above, with the 
full gravitational model employed. As with the RAAN=0 
case pair, most Fourier terms have been reduced from the 
point-mass gravity model’s case. 


Compression Curve Progression in RAAN at Summer Solstice, 
60 Degree Inclination, BN=100, Sweep in RAAN 
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Figure 6a. The compression curve changes in character with 
passage of the orbit near the subsolar point at the summer 
solstice. Here a 60° inclination orbit of 100 kg/m2 ballistic 

number was propagated every 15° in RAAN, ultimately 
lining up nearly through maximum density a few degrees 
east of the subsolar point when the initial RAAN at the start 
of the propagation was 120°. At RAAN = 300°(+180° of 
RAAN from 120°), the curve family reached a nearly 
mirrored amplitude profile. It illustrates the effect on 
symmetry of major perturbing effects in the atmosphere 
other than the equatorial bulge that causes the basic 
SIN(2©) curve. When all 22 runs are sequenced, the curves 
exhibit a simple harmonic “sloshing” with cycling RAAN. 
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Figure 6b. The variation in Fourier term amplitudes with 
RAAN for the sweep of a 60° orbit. Note that the 
20 amplitude term peaks well before the full compression 
curve does. 
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Figure 6c. The variation in Fourier term phases with RAAN 
for the sweep of a 60° orbit. 
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Figure 6d. When the first Fourier term in 1 is averaged at 
each value of © over all RAANS, a perfect sine wave 
emerges, representing the residue believed to be the 

J3 gravitational contribution. This slightly biases natural 
decays to occur in the northern (vastly more populated) 
hemisphere. 


From the Fourier amplitudes in Fig. 1 (especially in the 
90° and 270° RAAN cases) there is an apparent 
increasing trend in the peak 2© Fourier amplitude with 
increasing ballistic number. A special series of high 
ballistic number cases was run to explore this trend, with 
the peak amplitude at 90° inclination plotted in Fig. 7. 
From this we derive there is a simple quadratic trend in 
the symmetric compression effect as a function of 
ballistic number, while the asymmetric effects are 
increasingly less significant. 
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Figure 7. The trend of the first four Fourier term 
amplitudes for peak inclination orbits as ballistic 
number increases. The symmetric first and fourth terms 
show a simple quadratic growth, while the odd terms 
remain flat and comparatively small. 


It should be noted that it is only the symmetric terms that 
progress monotonically in phase, while the asymmetric 


terms present a less ordered structure (see Figs. 8 and 9). 


Trend of Fourier Cosine Phases vs. Ballsitic Number at 90° Inclination 
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Figure 8. The phases of the 20 and 40 Fourier (cosine) 
terms are the only monotonically changing curves in the 
significant components of the compression curves. 


The phases of the Fourier terms generally show a much 
noisier trend as a function of orbit inclination than do the 
amplitudes. A set of phases for one set of inclinations 


under common RAAN and BN conditions is shown in 
Fig. 9. As we have mentioned earlier, the compression 
curves do not exhibit particularly strong predictable 
structure in the lower inclinations. As in the high- 
ballistic-rnumber study, only the even-numbered 
(symmetric) terms show a monotonic and convergent 
behavior with increasing inclination. 
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Figure 9. The phases of the first four cosine terms of the 
compression curves for all ballistic numbers at 90° 
RAAN. There are four curves for each Fourier phase 
term: one for each ballistic number (50, 100, 150, and 
200 kg/m’). The curve for each successively higher 
ballistic number initiates farther lower on the graph 
than it did for the prior (smaller) ballistic number. Note 
the erratic behaviour of the asymmetric terms. (Phase 1 
and Phase 3, implying the offsets from ArgLat = 0° of 
the 10 & 30 cosine terms.) The jump appears 
coincident in inclination with the gap in the amplitude 
trend, evident in Fig. 4b). Only the 20 family is very 
well behaved over all inclinations. 
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Figure 10. The amplitudes of the first four Fourier 
terms of the compression curves, each averaged across 
four orthogonal RAANS. As with Fig. 9, each of four 
parallel curves starts lower than the curve for the next 
lower ballistic number. 


4. DISCUSSION 


The physics of decay biasing are difficult to discern 
empirically in highly nonlinear and _ interacting 
gravitational and atmospheric density domains that are 
each modelled with many orders of comparable-scale and 
often competing perturbations. Even the 30,000+ 
simulations of the current study are insufficient to isolate 
all the pure physics, although some key effects are 
emerging. 


The authors propose that the first four Fourier terms are 
sufficiently well-behaved in amplitude and phase that 
they can be used to generate a representative compression 
curve for any selected ballistic number and orbit 
inclination. The most problematic of these is the 
10 term, which as seen in Fig. 6d is likely harboring real 
J3 gravitational effects observed by Fremeaux, et al. [3] 
along with the diurnal atmosphere effects. While the 
diurnal effects are large, these should be averaged out 
when estimating expectation of casualty for decays 
sufficiently far in the future that specific dates and 
RAANSsS are completely unknown. 


BN= 50 ~~ BN= 100 
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Figure 11. The average of the first Fourier term over 
four orthogonal RAANSs at every value of ArgLat reveals 
a family of regular sine curves similar to that seen in 
Fig. 6d, whose amplitude and general regularity both 
grow with ballistic number. Each curve is a different 
inclination orbit, and all four plots are to the same 
scale. The curve amplitudes peak near 40° inclination 
and invert at the highest inclinations, but generally 
indicate a small preference (0.06 peak, vs 0.35 peak 
value in the two terms) to enter in the Northern 
Hemisphere for most orbits. 


Although only four RAAN’s were explored at the 
equinox (vs. the 22 at the solstice for one inclination 
shown in Figs. 6a-6d) the data can be similarly averaged 
to reveal (with some inherent noise) that the J3 effect 
seen in Fig. 6d is clearly at work. Fig. 11 is a compilation 
of the residue (average) curves of the 1@ term as a 
function of ArgLat at the equinox for each of four 


ballistic numbers 50, 100, 150, and 200, clockwise from 
upper left. 


The 30 phase term is sufficiently well behaved across all 
ballistic numbers and all inclinations below 85°, and its 
amplitude sufficiently small that it should be included, 
although its phase seems to be problematic. As 
evidenced by the point-mass_ gravitational study 
illustrated in Figs 5a-5d, the 3© term is expected to be 
atmospheric in nature. Fortunately, where the phase 
becomes noisy and suspect, the amplitude of the 
correction falls to near zero (perhaps an initial reason for 
the noise in estimating phase). 


4.1 The Model 


The authors continue to refine the model and the general 
understanding of the latitude biasing phenomenon in 
natural decays. The data presented has been used to build 
an approximate model of the compression curve for any 
set of conditions (see below). However, the cumulative 
errors in curve-fitting and then interpolating each Fourier 
term indicate that the most precise answer for any 
planned decay will be to perform a full simulation for the 
exact ballistic number and inclination of the vehicle in 
question, averaged over a continuum of RAANS (i.e., 
every 15°) on a date near the equinox. This does not take 
long using the analytical tools developed for this study. 


The next best model is to linearly interpolate the 
8 Fourier coefficients in phase and amplitude for any 
given ballistic number and inclination between the 
available 288 data sets. This model is what will actually 
be rolled into screening software and quick-study tools, 
as it is a tabular lookup that can run well in a spreadsheet. 
This method takes many thousands of data elements that 
cannot be presented easily in a publication, and which 
will be under ongoing update and refinement. 


The expedient, but least accurate approximate model 
presented here is an empirical derivation that—while far 
from perfect—is substantially more accurate than the 
prior assumption of a fully random decay along the 
ArgLat, and is useful for broad parametric studies to 
determine peak areas of concern in calculating E.. This 
model uses a set of “eyeball fit” empirical curves to 
model the change in coefficients. 


The model is built as follows: 
Each Fourier term Fy of Nth order has an Amplitude Ay 
and a phase offset ®y to create a term of the form: 


Fy=An*Cos (N*0+@y) 
In the equations that follow, BN is the ballistic number 


in kg/m’, and INC is the inclination of the orbit in 
degrees. 


a) F o=An= 1 


b) The first Fourier term represents the averaged 
first-term contributions over all RAANs for any 
given ballistic number. The resulting average is 
expected to be the effect of J3 gravitational 
perturbation as the principal 1© driver into the 
nonlinear atmosphere. Until further simulations 
reveal more detail, the evident simplest trend 
from the preliminary 4-RAAN average is a 
family of curves with the same phase: 

a. The amplitude of the first Fourier Term 
is estimated as the empirical fit: 
A =.0418+ 
.000125*BN*Sin((5.86*INC)+227 °) 
i. Note: the amplitude goes 
negative 
b. @&,=0° 
i. forlow BN 
@, = 100 *sin(2*INC+40°) 
ii. for BN >50 


c) The dominant second term F2 follows a sin? fit 
to inclination: 
a. A2=.0498 +0 .2576*sin3(INC) 
b. @2 = 192°*BN -170 
i. for INC >45° 
D2 = 70°*BN 
i. for INC <45° 


d) The F3 and Fy terms should most correctly be 
estimated from the averaged curves in all 
RAANs. This work has not been completed, 
and will be better done once a broader set of 
RAANs have been studied. The two terms are 
estimated from trend data at all four RAANs as 
follows: 

F 3. 
a. A3= .0334 -.0003*INC 
b. 3 =(BN)?.-60° 


e) The Fs term follows 
a. A4g=.0206-(8.5E-7)*(400-BN)*INC 
i. for INC < 55° 
A4=0.0206-(4.683E-5)*BN 
ii. for INC >55° 
b. Dy =130 - 
(BN/50)*(INC-50 °) 
i. for INC < 55° 
@4 = (135-BN/8)-INC 
i.  forINC > 55° 


0.78*BN + 


The fits of these respective parametric representations to 
the averaged Fourier terms is illustrated in Fig. 12. 
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Figure 12. The empirical fits described in this section to the 
Fourier coefficients for compression curves at each 
inclination, averaged across four orthogonal RAANS at the 
Vernal Equinox. The four curves in each plot represent the 
coefficients for the four studied ballistic numbers in 
ascending order left to right. The jumps in the 20 and 
30 phase curves do not lend themselves to simple models, 
and reconstructed curves in the regions of these breaks can 
show obvious qualitative differences. 


Two representative fits of the empirical curve to the 
averaged scenario and to representative original 
simulation data are shown in Figs. 13a and 13b for orbits 
not situated within 5° of a break point in the empirical 
curve sets. 
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Figure 13a. Empirical (dashed) fit to the 4-RAAN average 
(solid) model for moderate BN and inclination. Data trace 
(points) from 270° RAAN is included to demonstrate how 
actual conditions will typically vary from the average on any 
given day. 


Generic Model vs. 4-RAAN Avg for BN=200 and Inc =85 
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Figure 13b. The same comparison for high BN and 
inclination. The fit to the 4-RAAN average is excellent, but 
the statistical noise for any particular day of entry (in the 
form of the diurnal variation) is more evident. 


4.2 Final E- Calculation 


With the averaged concentration factor available at each 
latitude, the E, for a particular decay can finally be 
calculated with a slight modification to the traditional 
approach. The ArgLat at 90 km is incremented in small 
steps from 0° to 360°. The concentration factor is 
calculated for the ArgLat at 90 km for the specific orbit 
inclination and spacecraft ballistic number. A necessary 
pad of a few degrees is added to the ArgLat (variable with 
the ballistic number of the object in question) to 
propagate the debris farther downrange to the ground. 
(The calculation of this pad is beyond the scope of this 
paper.) The physical latitude on the planet associated 
with the adjusted impact ArgLat (Qimpact) is calculated by 
the simple equation: 


Latitude = ARCSIN(SIN(INC)*SIN(@impact)) (1) 


The population in any given latitude band is extracted 
from the gridded population of the world, with a model 
of population growth applied. In the critical adjustment, 
this population is divided by the concentration factor 
associated with the 90-km ArgLat © to give an adjusted 
risk. The binomial correction that must be applied to 
correct for a normalized Probability Function is 
P=(1/C())-(1-C(@))?. Note that the impact latitude will 
be encountered twice during this process (ascending and 
descending), with a different concentration factor 
applied. Each incremental risk is multiplied by the size 
of the incremental step in ArgLat divided by 360°, such 
that the sum of all steps is the total risk over the orbit. 


The resulting mean population per square kilometer can 
be mapped as a function of inclination and ballistic 
number, as shown in Fig. 14, illustrating the difference in 
resulting risk under the new model, compared to the prior 
randomized-entry assumption. Expectation of casualty is 


a linear multiple of the mean exposed population per 
square kilometer, scaling as well with Debris Casualty 
Area. The ratio of the average populations under various 
inclination orbits is a weak function of ballistic number, 
as shown in Fig. 15. 
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Figure 14. Weighted vs unweighted 2020 population under a 
90° orbit 
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Figure 15. Ratio of mean population density under 
compression curve vs. the assumption of uniform decay 
around the orbit, as a function of orbit inclination for four 
ballistic numbers. A Oinpact of 13.5 degrees downrange after 
the 90km altitude is included in this calculation. 


4.3 Comparison with Prior Models 


Prior work [1] identified a sinusoidal perturbation of the 
decay probability as a function of latitude, and implied a 
full symmetry of the effect. On an average case, this is 
approximately true, but for any particular orbit, the 
present work now shows the decay to be statistically 
clustered more strongly on the approach to the equator, 
as explained above. The physics of this has been 
discussed in [2], and there are significant ramifications 
for exceptionally low-dV intentional tactical de-orbit 
scenarios. In natural decays which usually occur over an 
indeterminate, broad timescale, many of these subtleties 


largely average out, but still have influence on the entry 
statistics. 


In the current work, the representation has been shifted to 
a function of Argument of latitude, and not latitude, to 
achieve a more universal formulation, and to search for 
underlying physics. The current work also recognizes 
and accommodates remaining asymmetries in statistical 
decays, and refines the model for ballistic number 
differences. 


The prior method of estimating E,. incorporated a 
rudimentary concentration factor based upon how long a 
circular orbit dwells over any given latitude. Typically, 
this correction is applied symmetrically around the 
equator in an algorithm to adjust the latitude-averaged 
population bands in strips only a few tens of kilometres 
wide to account for this circular “dwell time”. In the 
method described within this paper, this correction is not 
applied independently, substituting instead a localized 
concentration factor in tiny numeric integration steps in 
ArgLat, applied and integrated band by band. Note that 
the prior method actually places stronger likelihood of 
decay at the latitude extremes, and minimizes the effect 
at the equator, in complete opposition to the effect 
explored here. The effective weighting of population is 
evident in Fig. 14. 


To conclude the forward work described in the authors 
prior study [1], the influence of F(10.7) solar flux 
variations was shown in [2] to be negligible in the critical 
final eight orbits when the decay rate is geometrically 
increasing. This is because the expansion of the 
atmosphere results from radiant energy deposited above 
the altitudes where the lowest density scale heights (and 
thus, maximum effect of perturbations) occur. 


Conversely, the global geomagnetic index (Kp) was 
shown to be the driving factor to disturb low altitude 
atmospheric density. Geomagnetic index (updated every 
3 hours) determines the energy deposition of charged 
particles in the altitudes driving the final few orbits of 
decay. From a statistical standpoint, the Kp variation is 
rapid and random enough to only add noise to the more 
consistent geodetic effects of the wall of air and of 
J2 perturbations. That said, the effect of Kp is significant 
when attempting a drag-centric targeted decay on any 
particular day. 


5. CONCLUSIONS 


The statistical clustering of natural decays towards the 
equator has been reduced to an approximate parametric 
representation accounting for ballistic number, 
inclination, and gravitational perturbations, useful for 
parametric studies and potentially for insight into the 
physics. Some of the corrections to prior practice can 
lead to localized E, changes of up to 36% over the 


random-decay assumption, and net integrated E,’s that 
are 18% or more different than prior practice would 
estimate. There is strong evidence that a slight bias of 
natural decays occurs in the northern hemisphere due to 
the J3 gravitational effect. 


However, the chaotic nature of some parameters 
indicates that best precision will come from a dedicated 
simulation of each decay’s ballistic number and orbit 
inclination (and especially, if known, the date and range 
of likely RAANSs for the expected decay), rather than the 
inherently “fuzzy” approximation to the specific case in 
a set of global parameters that seek to optimally fit a 
broader range of cases. 


The clustering is most conveniently rendered as a 
function of the ArgLat, which can then be transformed to 
latitude. With the local compression factor available, it 
is a minor adjustment to existing look-up tools to 
generate a statistically- weighted E, that accounts not only 
for the inhomogeneity in the population density with 
latitude, but as well for the inhomogeneity in the 
distribution of decays around the orbit. 
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